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Abstract 

Multi-configuration range-separated density-functional theory is extended to the time-dependent 
regime. An exact variational formulation is derived. The approximation, which consists in combin- 
ing a long-range Multi- Configuration- Self- Consistent Field (MCSCF) treatment with an adiabatic 
short-range density-functional (DFT) description, is then considered. The resulting time-dependent 
multi-configuration short-range DFT (TD-MC-srDFT) model is applied to the calculation of singlet 
excitation energies in H2, Be and ferrocene, considering both short-range local density (srLDA) and 
generalized gradient (srGGA) approximations. In contrast to regular TD-DFT, TD-MC-srDFT can 
describe double excitations. As expected, when modeling long-range interactions with the MCSCF 
model instead of the adiabatic Buijse-Baerends density-matrix functional as recently proposed by 
Pernal [K. Pernal, J. Chem. Phys. 136, 184105 (2012)], the description of both the l 1 !? doubly- 
excited state in Be and the state in the stretched H2 molecule are improved, although the 
latter is still significantly underestimated. Exploratory TD-MC-srDFT/GGA calculations for fer- 
rocene yield in general excitation energies at least as good as TD-DFT/CAM-B3LYP, and superior 
to wave-function (TD-MCSCF, symmetry adapted cluster-configuration interaction) and TD-DFT 
results based on LDA, GGA, and hybrid functionals. 
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I. INTRODUCTION 



Time-dependent density-functional theory (TD-DFT) (TJ |2] is used routinely nowadays 
for computing electronic excitation energies and transition properties of molecules and solids. 
Even though low-lying excitation energies often can be obtained sufficiently accurate at a 
relatively low computational cost, adiabatic TD-DFT using pure functionals usually fails in 
describing, for example, charge transfers and double excitations. While charge transfers often 
can be modeled adequately with range-separated hybrid functionals [3H7], where the long- 
/short-range decomposition of the electron-electron repulsion is used only for the exchange 
energy, the double excitations remain problematic for all standard exchange-correlation func- 
tionals as long as the adiabatic approximation is used for the time-dependent exchange 
and correlation density-functional [8HTT]. This statement holds for hybrid functionals even 
though non-adiabatic effects can be taken into account through the exact exchange func- 
tional, due to the fact that the Kohn-Sham (KS) orbitals are non-local (in time) function- 
als of the density [T^HH] . The accurate description of charge-transfer states in molecular 
complexes is of particular interest because of their prominent role in for example organic 
electronics [15] and dye-sensitized solar cell applications [IE]. To this end, more sophisti- 
cated functionals such as LRC-BOP [Lj], CAM-B3LYP [18J, as well as LC-w-PBE [H] 
and other types of range-separated hybrids have been developed [201422] . 

In order to further improve standard TD-DFT, in particular for the description of double 
excitations, Pernal [23J recently proposed to combine it with time-dependent density-matrix 
functional theory (TD-DMFT) by means of range separation. While describing the long- 
range part of the electron-electron interaction with the Buijse-Baerends (BB) density-matrix 
functional, the short-range interaction is treated within the local density approximation 
(srLDA), both within the adiabatic approximation. The author has tested on H 2 and Be, 
and she found that the combined method performs much better than standard TD-LDA and 
TD-DMFT-BB for many excitations, provided the range separation parameter is properly 
chosen. However, large errors were obtained for the l^D double excitation in Be as well 
as for the l 1 ^ state in the stretched H2 molecule. As an alternative to the long-range 
adiabatic TD-DMFT-BB treatment, we propose in this work to describe the long-range 
interaction at the time-dependent Multi- Configuration Self-Consistent Field (TD-MCSCF) 
level, extending thus the MC-srDFT [21] model to the time-dependent regime. 
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The paper is organized as follows: after a short introduction to time-independent range- 



separated DFT (Sec. II A), an exact and variational time-dependent formulation is presented 



in Sec. II B The approximate linear response TD-MC-srDFT scheme, which consists in 
describing the long-range interaction at the TD-MCSCF level and the short-range interaction 



within adiabatic TD-DFT, is then introduced and discussed in Sees. |II C| and |IID[ The new 
scheme is applied to two widely varied types of atomic and molecular systems. We first study 
the paradigm systems H 2 and Be to illustrate fundamental benefits of the TD-MC-srDFT 
approach by comparing to other suggestions for how to go beyond Kohn-Sham DFT. We 
further demonstrate the performance of the new TD-MC-srDFT scheme in an investigation 
of the low-lying singlet excitations of ferrocene. Thus, following the computational details 



(Sec. Ill), numerical results obtained for singlet excitations in H2 (Sees. IV A and IV B), Be 



(Sec. IV C) are discussed. In Sec. IV D we elaborate on valence- and charge-transfer singlet 



excitations in the transition-metal compound ferrocene before drawing conclusions in Sec.fVj 



II. THEORY 



This section deals with the rigorous formulation of variational range-separated TD-DFT 
models where the range separation is used for both exchange and correlation energies. Intro- 



ducing first range-separated DFT in Sec. |II A[ its extension to the time-dependent regime is 



then presented in Sec. |II B In particular, the derivation of an exact variational formulation 
is detailed, in the light of Vignale's recent work [251 EH] on the Runge and Gross varia- 
tional principle in TD-DFT. For practical calculations, an adiabatic variational formulation 
is also given. The approximate TD-MC-srDFT model, where the long-range interaction is 



described within the TD-MCSCF approach, is then derived in Sec. II C using Floquet theory. 



The interpretation of the TD-MC-srDFT poles as excitation energies is then discussed in 

sec. una 

The theory in this section is completely general in the sense that it is valid for any reasonable 
range-separation of the two-electron repulsion, however, we will for simplicity describe the 
theory with the most commonly used choice for the range-separation: 



^e r /(r 12 ) = erf(/iri 2 )/r 12 , 



wff{r 12 ) = erfc (/ir 12 )/r 12 = l/r 12 - w l *f{r 12 ), 



(1) 
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which depends on the parameter /i. We note that with this choice one has a continuous 
class of different range-separations when varying \x from zero to infinity: in the fi = limit 
the long-range part is zero, w^°(r 12 ) = 0, while in the \x = +00 limit the short-range part 
is zero, ?«^ +00 (7*12) = 0. 

A. Range-separated density-functional theory 

In multi-determinant range-separated DFT, which in the remainder of this paper is simply 
referred to as short-range DFT (srDFT), the exact ground-state energy of an electronic 
system can be expressed as 

E = mm {<*|T + Wl r f + V ae m + , (2) 

where T is the kinetic energy operator, denotes the long-range two-electron interac- 

tion, V ne is the nuclear potential operator, and ^HxcN denotes the //-dependent short-range 
Hartree-exchange-correlation (srHxc) density-functional which describes the complementary 
short-range energy f27\ f2E\ . The minimizing wave function in Eq. Q fulfills the following 
self-consistent equation: 

(f + W 1 :/ + V ne + V££[nw]j = 

(3) 

where n(r) is the density operator. It is readily seen from Eqs. @ and Q that srDFT 
reduces to standard Kohn-Sham (KS) DFT and pure wave function theory (WFT) in the 
fi = and /1 — » +00 limits, respectively. A combined WFT-srDFT approach is then obtained 
for intermediate fi values. The approximate Hartree-Fock-srDFT (HF-srDFT) scheme con- 
sists in restricting the minimization in Eq. ^ over single-determinant wave functions. It 
is very similar to range-separated hybrid schemes where the range separation is only used 
for the exchange energy [3H7]. The two approaches differ by the fact that HF-srDFT does 
not describe long-range correlation effects since, in srDFT, the range separation is also used 
for the correlation energy. In the case of multi-configurational systems, the minimization 
in Eq. ^ should be performed over MCSCF-type wave functions instead. This scheme 
is referred to as MC-srDFT [2^1 129] in the following. Note that correlated methods such 
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as second-order M0ller-Plesset perturbation theory (MP2) [30] , Random-Phase Approxima- 
tion (RPA) j3U [32], Coupled-Cluster (CC) [33] and second-order n-electron valence state 
perturbation theory (NEVPT2) [31] have also been applied in this context, using the non- 
variational formulation of srDFT which is given in Eq. ([3]). It is important to mention 
that, when describing the long-range interaction at the post-HF level and the short-range 
interaction with a short-range density-functional, there is no risk of double counting electron 
correlation effects, which is essential for a combined post-HF-DFT approach. In this respect, 
using a range-separation of the two-electron repulsion is very appealing. On the other hand, 
the performance of post-HF-srDFT models will depend on how correlation effects can be 
split into long- and short-range contributions. This is clear in van der Waals systems in 
which, for example, MP2-srDFT or CC-srDFT models perform relatively well [301 [33]- In 
multi-configurational systems, however, it is in general not possible to interpret static and 
dynamical correlations as long- and short-range ones, respectively. As a result, even though 
MC-srDFT performs better than regular DFT in stretched molecules for example [24], the 
approximate short-range density-functional part of the energy is usually not accurate enough 
since it has somehow to describe a part of static correlation [35]. In this respect, approximate 



short-range functionals that have been developed so far (see Sec. Ill) need improvements 
and work is currently in progress in this direction. 

In this paper, the time-dependent regime will be explored in details only for the variational 
HF-srDFT and MC-srDFT methods. The extension to non-variational srDFT schemes, 



which is briefly mentioned in Sec. II B is left for future work. 



B. Variational principle in the time-dependent regime 

The extension of exact srDFT to the time-dependent regime could in principle be achieved 
when considering the time-dependent Schrodinger equation of an auxiliary long-range inter- 
acting system [23], according to the first Runge and Gross theorem PQ, that is without using 
a variational principle. If the long-range interaction is treated approximately, at the MC- 
SCF level for example, the Ehrenfest theorem [3E] could then be applied in order to obtain 
the time-dependent MC-srDFT response functions. As an alternative, we explore in this 
section the possibility of formulating exact time-dependent srDFT variationally, using the 
quasienergy formalism [37] . As pointed out in Ref. [3H], the latter is arguably more attractive 
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than the Ehrenfest method in that it provides a unified framework for applying variational 
and non- variational long-range post-HF methods, by analogy with time-independent theory, 
to which it naturally reduces in the limit of a static local potential As an additional 
advantage, the permutational symmetries (with respect to the exchange of perturbation op- 
erators) is manifest in the quasienergy method. 

In order to obtain a time-dependent extension of Eq. ([2]), we use in the following the recent 
work of Vignale [25] [26] . Let us consider the action integral p] expression 

fi[*]= F Q[*](t) dt, 

(4) 



(*(t)\T+w„+v(t)-a\*(t)) 



Qim) 



(*(*)!*(*)> 

which is defined for a given time-dependent wave function ^(t) and a given time-dependent 
local potential V(t) = J dr v(r,t)h(r). The operator W ee denotes the regular two-electron 
repulsion. Note that, in Eq. Q, ^(t) is not assumed to be normalized so that, when con- 
sidering infinitesimal variations of the time-dependent quasienergy Q in the following, 
no normalization constrain will be needed. For time-independent local potentials and wave 
functions, the time-dependent quasienergy reduces to the usual energy. In Floquet theory, 
which is considered in Sec. |II C the action integral over a period T equals what is usually 



referred to as a quasienergy [37J EH] multiplied by T. In the general (non-periodic) case, the 
quasienergy would be equal to the action integral divided by (ti — 1 ), which is a constant. 
For simplicity, we will refer to Q[^} as quasienergy even though it is an action and not an 
energy, strictly speaking. The time-dependent Schrodinger equation 

f + W ee + V(t) - d) Mt)) = Q[¥](f)|*(f)>, (5) 



dt / 

where \l/(io) is assumed to be normalized, is then equivalent to the variational principle 
based on the quasienergy 

6Qm = -i\(^(t)\5^(t)}} t \ (6) 

L J to 

where variations \l/(t) — > \l/(t) +S^(t) around the exact solution *&(t) are considered. Using 
the boundary conditions 5^(t ) = = 0, the variational principle can be simply written 

as 

6Q[V] = 0. (7) 
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This stationary condition is convenient for deriving approximate time-dependent response 
properties based on variational methods such as HF and MCSCF [37] . The time-dependent 
extension of both HF-srDFT and MC-srDFT schemes can then be achieved along the same 
lines, provided that time-dependent srDFT can be expressed in terms of a variational prin- 
ciple, which is actually not trivial [25, 26J. This point is addressed in the rest of this section. 
According to the first Runge and Gross theorem pQ, for a fixed initial wave function *&(to), 
the following density-functional quasienergy can be defined: 

Q[n} = B[n} + jf dt J dr v(r,t)n(r,t), (8) 

where the universal functional 

B[n] = jf 1 dt (*[n](t)\f + W ee - i^\*[n)(t)), (9) 

is calculated with the time-dependent wave function \l/[n](t) associated to a fully interacting 
system whose time-dependent density equals n(r,t): 

f + W ee + V[n)(t) - i^j Mn]{t)} = Q[n](t)\*[n)(t)}, 

(10) 

V[n](t) = I dr v[n](r, t) h(r). 



The density n(r,t) associated to the exact solution of the time-dependent Schrodinger 
Eq. ^ can then be obtained when applying the following variational principle 

5Q[n] = -i<%](t0l^[n](*i)>, (11) 



where the boundary condition <5\I/ [fi] (to) = was used. As pointed out by Vignale 
though the boundary conditions 5n(r, to) = 5n(r,ti) = are fulfilled, the condition 5Q[n] = 
is not. Indeed, the wave function ty[n + 5n](t) = ty[h](t) + Sty[h](t) accumulates the change 
of the density 5n(r, t) over the time interval [to ; ^i] : 

f* 1 f Sty 
<ftt[n](fi) = J dt J dr [w](ti)<fa(r,t) 

+ o|° (12) 



so that the variational principle in Eq. (jll[) can be rewritten as 

SQ 



As in the time-independent regime (see Sec. II A), the universal functional B[n] defined in 
Eq. ^ can be split into long-range and short-range parts: 



n\ 



(14) 



to 



dt(^[n](t)\f + W 1 :/ - ij^>](t)> 5 



where ^ M [n](t) is the time-dependent wave function associated to a long-range interacting 
system whose density equals n(r,t): 



f + ti£* + V>](t)-i^) !*>](*)> 



Q»[n](t)\*»[n](t)), 



(15) 



V"[n](t) = / drv'*[n](r,t) n(r). 

Note that the initial state \l/ M [ra](to) is fixed and equal to the auxiliary long-range interacting 
wave function \l/ M (to) whose density equals the one of the real fully interacting initial wave 



function ^/(t ). The universal long-range functional defined in Eq. (14) can be expressed in 



terms of the auxiliary long-range interacting time-dependent quasienergy as follows 

B ir )M [ n ] = f 1 dt (Q^[n}(t) 
J to V 



— J dr i7 M [n](r,t) n(r,i)j. 



(16) 



so that, using the simplified expression 



dt 5Q"[n](t) = -i(V»[n](t 1 )\5W i [n]{t 1 )) 



(17) 



- J dt I dr 5v fl [n)(r,t)n(r,t), 



as well as Eqs. <\8h and (14), the variational principle in Eq. (13) provides the exact expres- 



sion for the local potential which reproduces the time-dependent density n(r, t) of the fully 



interacting wave function from a long-range interacting one: 

v"[n){r,t) = v(r,t) + 



Hxc 




5n(r, t) 



+i (* |fil(ii) *^)' fi|(ti) >- (18) 

It is important to keep in mind that, as further discussed in Sec. |IID the time evolution 
of the auxiliary long-range interacting system yields the time-dependent response of the 
exact density but not the response of the exact wave function. In the particular case of a 
periodic perturbation, which is considered in the rest of this work, both fully and long-range 
interacting wave functions remain unchanged after a period T: 

*[n](to + T)=#[n](t )=*(*o), 

(19) 



so that, for t\ = to + T, the last two terms in the right-hand side of Eq. (18) are equal 



to zero. In addition, if the adiabatic approximation is used for describing the short-range 
interaction, the srHxc density-functional quasienergy expression is simplified as follows: 



Q 



Hxc 



dtE s ^ c [n(t)i 



(20) 



to 



where E^^.[n] is the time-independent srHxc density-functional introduced in Eq. ([2]), and 



the time-dependent potential in Eq. ( 18 ) becomes 



5n(r) 



(21) 



Combining Eqs. ([3]), ( 15 ) and (21 ), we conclude that the time-dependent density n(r, t) of the 
real fully interacting system can be approximated from an auxiliary long-range interacting 
one whose wave function \I/ M (t) fulfills, within the short-range adiabatic approximation, 

.8 



T + Wlr + V(t)+V^[n Mt) ] 
Q»{t)\W(t)), 



dt 



(22) 
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which is equivalent to the stationary condition 

m] = 0, (23) 

where the wave-function-dependent range-separated quasienergy is defined as 

"" ^(t)\f + W^ + V(t)-i§- t \^(t)) 



= r dt 

Jtn 



mmt)) 

(24) 



+ j dt E s ^ c [n m }. 

J t 

Note that the time-dependent KS equation, as formulated within the adiabatic approxi- 



mation, is recovered from Eq. (22) when /i = 0, while the fi — > +oo limit corresponds to 
the time- dependent Schrodinger equation. When < /i < +oo, a rigorous combination of 
density-functional and wave function theories is obtained in the time-dependent regime. As 
shown in the following, a multi-configuration extension of regular TD-DFT can then be for- 
mulated in this context when describing the long-range interaction at the MCSCF level. Let 
us mention that non- variational methods such as CC could also be merged with TD-DFT, 



using Eq. (22) in combination with a Lagrangian formalism |37j. 



C. Multi-configuration range-separated TD-DFT based on Floquet theory 

We work in this section in the framework of Floquet theory [37] where the time- dependent 
perturbation is periodic: 

N 

x k=-N 

= -jr, (25) 
V x = dr v x (r)n(r), 



and the quasienergy in Eq. (24) is calculated over a period f* 1 dt — > Jo dt. Since the long- 



range interaction in Eq. (22) is treated explicitly, the exact time-dependent wave function 
^(i) is a multi-determinant one. As an approximation, we consider the following MCSCF- 
type parametrization consisting of exponential unitary transformations 



\W(t)) = e^'V^TO, (26) 
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where ^/q denotes the unperturbed time-independent MC-srDFT wave function and 



K 



(*) = E + e" 



(27) 



The singlet excitation and state-transfer operators are defined as follows 

<?1 = 4*? = «LAa + 4 A/3! p > q, 



(28) 



R] 



Note that the TD-HF-srDFT scheme is a particular case of Eq. (26), where the unperturbed 



MC-srDFT wave function would be replaced by the HF-srDFT determinant, and orbital 
rotations only would be considered. The resulting linear response equations would then be 
formally identical to standard TD-DFT equations based on hybrid density-functionals 



Returning to the multi-configuration case, the TD-MC-srDFT wave function in Eq. (26) is 
fully determined by the Fourier component vectors 



A(w,) 



Ki(-Wl) 



(29) 



for which we consider in the following the Taylor expansion through first order: 



N 



k=-N,x 



de x (uj k ) 



+ 



(30) 



(31) 



Rewriting the variational condition in Eq. (23) as follows 

the linear response equations are simply obtained by differentiation with respect to the 
perturbation strength e x [39| 130]: 

d dQ" 



de x (uj k ) 9At(-wj) 



0. 



(32) 
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According to Eq. (24), the quasienergy can be decomposed as follows: 



Hxc> 



(33) 



where the purely long-range MCSCF part equals 



Q 1 ^ = / dt (tf"(t) | f + Wl v f + V(t) I 



+ 1 dt(*"(t) \ -i^l*"(t)>, 

and the purely short-range DFT contribution is written as 



(34) 







(^(t)|n(r)|^(t)}. 



(35) 



The terms arising from the derivatives of the former can be computed with a regular linear 
response MCSCF code [36, 41 j , using long-range two-electron integrals, and thus do not re- 
quire additional implementation efforts. On the other hand, standard TD-DFT codes cannot 
be used straightforwardly for computing the srDFT terms since the density is now obtained 
from a MCSCF-type wave function instead of a single KS determinant, terms describing 
how the density changes when the configuration coefficients change are also needed. The 
srDFT contributions to the linear response equations can be decomposed as follows 



d 



dQ s 



de x (u k ) dA^(-ui) 



d 



a 




{^(t)\V^[n^(t)} dt 



+—^2S(Lo m + Lo n ) J drdr' K\ 



(36) 



Hxc V- 1 



r, r 



x A^-co m )n^(v)n^(r')A(co n ) 
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where if^(r, r') 



5 2 E^J8n{v)8n{v')[n^ 



denotes the srHxc kernel calculated for the 



unperturbed density and the gradient density vector equals 



n 



[iW 



(*»,™(r)» 
<*SI[$Uto]TO 



(37) 



The linear response Eq. (|32|) can thus be rewritten as [36l 



de x (u k ) 



(38) 



where the Hessian is split as follows: 



E [2]» = E [2]» + K s*£ 



(39) 



The MCSCF-type Hessian is based on the auxiliary long-range interacting Hamiltonian 
f + + V ne + V££[n 9 »], that is used as H Q in Eqs. (9) and (10) of Ref. jH], and the 
srHxc kernel contribution is defined as 



-"■Hxc 



drdr' K^(r, r' )// 



Mm, 



rm 



M/4/V 



(40) 



Both E [ q ]i1 and S [2] ^ matrices (see Eqs. (7) and (8) in Ref. [H]) are built from the un- 
perturbed MC-srDFT wave function \1/q and the gradient property vector Vx 1 ^ equals 
J dr v x (r)n^(r). Excitation energies can thus be calculated at the MC-srDFT level when 
solving iteratively 



E Wp_ uS 12]A X 



(41) 



which can be considered a multi-configuration extension of the Casida equations [2] . 



D. Interpretation of the TD-MC-srDFT poles 



The linear response Eq. (41 ) computes the poles of the TD-MC-srDFT wave function. As 



already mentioned in Sec. |II B[ a solution can only be interpreted as an excitation energy if 
it is a pole of the TD-MC-srDFT density. In this respect, the double excitation recovered in 
the n = limit of a simple two-state TD-MC-srDFT model (see Appendix) is not physical. 
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It is in principle less problematic for non-zero fi values, even small ones, since the ground- 
state MC-srDFT wave function becomes multi-determinantal. This is analyzed further in 



Sees. IV B and IV C In order to obtain a smoother connection to regular TD-DFT when 
/j, — > 0, an effective orbital rotation vector 



ic(ui) 



(42) 



such that 



n[l]/4l 



n 



nW^(r) 



(43) 



could be introduced and the TD-MC-srDFT linear response Eq. (38) reformulated in terms 



of dK{— LJi) / dsJcUk) ■ Work is currently in progress in this direction. 

o 



III. COMPUTATIONAL DETAILS 



The TD-MC-srDFT linear response Eq. (38) has been implemented in a development 



version of the DALT0N2011 program package [42J. Calculations have been performed with 
spin-independent short-range functionals, considering both local density (srLDA) and gener- 
alized gradient (srGGA) approximations. In the former case, we used the srLDA functional 
of Toulouse, Savin and Flad [31 H3]. We used, as srGGA functional, the Perdew-Burke- 
Ernzerhof-type (PBE) functional of Goll, Werner and Stoll [33], which is denoted srPBEgws 
in the following. Basis sets are aug-cc-pVQZ [H] for both H2 and Be systems. We fur- 
thermore performed all-electron TD-MC-srDFT linear response calculations of ferrocene 
Fe(CsH5)2 (iron frzs-cyclopentadienyl, FeCp2) examining excitation energies and oscillator 
strengths of the lowest singlet excited d-d and ligand-to-metal charge transfer states. The 
geometrical parameters for the eclipsed FeCp2 conformer (D^ symmetry), compiled in Table 
SI in the Supplementary Material |45j . were taken from the recent work by Coriani et al. 
[46] who carried out highly accurate geometry optimizations at the CCSD(T) level yielding 
close agreement with experiment and other available ab initio data [4TH50] . We employed 
triple-^ cc-pVTZ-DK basis sets for all elements [5U |52] where scalar-relativistic effects were 
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taken into account by means of the Douglas- Kroll-Hess second-order (DKH2) Hamiltonian. 
We did not account for spin-orbit effects which were found by Scuppa and co-workers [30] 
to be of only minor importance for the singlet excitation spectrum of ferrocene. The excited 
state manifold of FeCp2 was computed with the regular TD-HF and TD-MCSCF approaches 
as well as the combined TD-HF-srLDA, TD-HF-srPBEgws and TD-MC-srPBEgws models. 
The initial CASSCF optimization step was performed for both MC-srDFT and regular M C- 
SCF using the well-established CAS(10,10) active space with 10 electrons in 10 orbitals 
[IB] ) comprising the Cp-ligand 7r orbitals in addition to the Fe 3d4s shells. For analy- 
sis purposes, we carried out standard scalar-relativistic DKH2 TD-DFT calculations using 
pure LDA [53] and PBE [31] functionals, as well as the hybrid Becke three-parameter Lee- 
Yang-Parr B3LYP [55] and the Coulomb attenuated method CAM-B3LYP [18J functionals, 
within the adiabatic approximation. The /i parameter was set in all TD-HF-srDFT and 
TD-MC-srDFT calculations to 0.4 a.u. unless otherwise specified. This value relates to a 
prescription given in Refs. [2"4"1 135"] for an optimal treatment of short-range electron corre- 
lation in ground-state MC-srDFT calculations and should, in principle, be re-considered in 



the time- dependent regime. As illustrated in Sec. IV C the choice of // in TD-MC-srDFT 
calculations is important since it affects excitation energies significantly. This will be further 
analyzed in future work. 



IV. RESULTS AND DISCUSSION 

A. 1 S+ excited states of H2 

Our TD-MC-srDFT calculations have been performed for H 2 with the minimal la g la u 
active space which however can recover excitation energies very close to the FCI ones 
already at the TD-MCSCF level (see Fig. [I] (a)). Potential curves obtained with the srLDA 
and srPBEgws functionals were found to be on top of each other. Therefore only the formers 
are discussed in the following. As shown in Fig.[l](b), the TD-HF-srLDA and TD-MC-srLDA 
methods give, near the equilibrium H-H distance (about 1.4 a.u.), the same excitation ener- 
gies for the first four x £+ states when \x is set to 0.4 a.u. This is due to the fact that (i) the 
ground-state HF-srLDA and MC-srLDA wave functions are, in this case, almost identical [21] 
(ii) these excitations all correspond, predominantly, to single excitations from the la g to a u 



16 



orbitals [57J. The latters are indeed well described by orbital rotations, like in standard TD- 
DFT, with no need for a long-range multi-configuration treatment. Note that the excitation 
energies are, in this case, more accurate than the pure TD-LDA ones: using the FCI curves 
as reference, the error is approximately divided by two. Upon bond stretching, the difference 
between TD-HF-srLDA and TD-MC-srLDA excitation energies becomes significant for the 
4 1 X+ state already around 2.5 a.u. Indeed, while the former is associated to the singly ex- 
cited configuration (1<t 9 ) 1 (4ct u ) 1 , the latter corresponds to the doubly excited configuration 
(l(j u ) 1 (2(T 5 ) 1 . This double excitation is obtained as a single excitation applied to the doubly 
excited configuration (1<7 U ) 2 whose weight in the ground-state MC-srLDA wave function is 
less than 1% but not strictly equal to zero. For the same bond distance, the 3 X E^ state, 
which is dominated by the singly excited configuration (la 5 ) 1 (3cr M ) 1 , acquires also, but to a 
less pronounced extent than for 4 X E+, a doubly excited character [{la g ) 2 — »■ (lcr u ) 1 (2<j g ) 1 ]. 
The and 2 1 £+ states are mainly combinations of the singly excited configurations 

(lcr 9 ) 1 (lcr u ) 1 and (lcr 9 ) 1 (2a n ) 1 , which explains the smaller difference, compared to the 4 1 £+ 
state, between the TD-HF-srLDA and TD-MC-srLDA excitation energies. Around 5 a.u., 
the 4 1 E+ state obtained at the TD-MC-srLDA level is dominated by the doubly excited 
configuration (lcr u ) 1 (3o" 9 ) 1 while, at the regular TD-MCSCF level, both (la u ) 1 (3a g ) 1 and 
(la u ) l (2a g ) 1 configurations are important. This difference could be justified by the fact that 
long-range interactions only are described within the MCSCF model while short-range inter- 
actions are assigned to DFT. In other words, the ground-state wave function and its linear 
response are not expected to be the same at the MC-srDFT and regular MCSCF levels. 
As shown in Sec. |IIB this should be expected only for the densities. The large underesti- 
mation of the 4 X E+ excitation energy, by about 0.1 a.u., is mainly due to the approximate 
(adiabatic srLDA) potential and kernel used. It is important to notice that the explicit 
treatment, at the MCSCF level, of the long-range interaction enables a multi-configuration 
description, within a TD-DFT framework, of the excited states upon bond stretching. This 
is illustrated by the increasing doubly excited character [(lcx g ) 2 — > (la u ) 1 (2a g ) 1 ] of the 3 X S+ 
state when enlarging the bond distance from 2.5 to 5 a.u. This explains the increasing dif- 
ference between the TD-HF-srLDA and TD-MC-srLDA excitation energies. For the latter, 
the doubly excited character is induced by the single excitation [(lcr n ) 2 —> (lcx tt ) 1 (2cr 9 ) 1 ] 
applied, in the long-range MCSCF linear response calculation, to the configuration (lcr u ) 2 
whose weight increases in the ground-state MC-srLDA wave function (0.6% at 2.5 a.u. and 
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23% at 5 a.u.). Let us mention that, even though TD-HF-srLDA and TD-LDA methods 
cannot describe double excitations, they provide at 5 a.u. a more accurate 3 1 E+ excita- 
tion energy than the TD-MC-srLDA approach, an example of fortuitous error cancellation. 
Upon further stretching, both TD-LDA and TD-HF-srLDA potential curves indeed devi- 
ate increasingly from the FCI one, while the TD-MC-srLDA curve remains parallel to the 
latter. Such a fortuitous error cancellation does not happen for the 2*E+ and states, 
and the TD-MC-srLDA performs better than TD-LDA and TD-HF-srLDA methods upon 
bond stretching. At 5 a.u., the 2 1 £+ state is dominated by both singly excited (lo" g ) 1 (2cr u ) 1 
and doubly excited (la u ) 1 (2a g ) 1 configurations for which the MC-srLDA linear response 
coefficients are equal to 0.33 and 0.13 in absolute value, respectively. On the other hand, 
the l^Ej state is dominated by the singly excited (lcr 9 ) 1 (lcr u ) 1 ionic configuration which is 
included in the active space and thus described, not by orbital rotations, but by configu- 
ration rotations instead. As expected [57], the TD-HF-srLDA 1 X S^ potential curve does 
not exhibit a minimum, exactly as with TD-LDA, simply because the la g and la u orbitals 
become degenerate in the dissociation limit. On the contrary, the TD-MC-srLDA curve has 
a minimum. This clearly shows that describing the long-range part of the electron-electron 
interaction at the MCSCF level in a range-separated TD-DFT framework can improve re- 
markably the exchange-correlation kernel used in standard TD-DFT calculations. Note also 
that the TD-DMFT-srLDA (/i = 0.7) 1*£+ excitation energy of Pernal [23] (3.02 eV), which 
was computed for a bond distance of 10.0 a.u., is largely improved at the TD-MC-srLDA 
(fi = 0.7) level (8.45 eV) when compared to FCI (9.90 eV). All these results were obtained 
with the aug-cc-pVTZ [44J basis set. Even though we use a different srLDA functional, we 
can reasonably conclude that the large errors obtained at the TD-DMFT-srLDA level are 
essentially due to the long-range adiabatic BB approximation. Returning to the TD-MC- 
srLDA (fi = 0.4) results, the 1 X S^ excitation energy is still significantly underestimated for 
large bond distances. The srLDA potential and kernel should then be improved in order to 
obtain more accurate results. One can also notice that the FCI l 1 ^ potential curve min- 
imum (around 4 a.u.) is shifted to about 5.5 a.u. at the TD-MC-srLDA level. Indeed, up 
to bond distances of 4 a.u., there are no significant differences between TD-HF-srLDA and 
TD-MC-srLDA excitation energies. This is due to the relatively small fi = 0.4 value, which 
ensures that, at the equilibrium distance, most of the electron correlation is assigned to the 
short-range interaction and thus treated in DFT. The assignment of static correlation to the 
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long-range interaction appears clearly in the dissociation limit, that is when both la g and 
la u natural orbitals are singly occupied in the ground-state MC-srLDA wave function [24]. 
For intermediate bond distances, the range-separation of static and dynamic correlations 
is not clear anymore [35]. It means that the short-range potential and kernel should be 
accurate enough to obtain reliable excitation energies for distances where a part of static 
correlation is inevitably assigned to the short-range interaction. The situation is of course 
less critical for larger \x values but, in this case, more electron correlation must be described 
by the MCSCF which is not appealing, in terms of computational cost, for larger scale cal- 
culations. Using a multi-configuration short-range exact exchange energy expression [SSHBT)] 
while keeping \x = 0.4 a.u. is a possible alternative currently under investigation. Such a 
scheme would largely reduce self-interaction errors [21] in the ground-state MC-srLDA en- 
ergy and is expected to affect excitation energies through the improved short-range potential 
and kernel. 



B. excited states of H2 

The TD-MC-srLDA method was also applied to the calculation of the first four 
states of H 2 , setting \x to 0.4 a.u. and using the minimal \a g \a u active space. The latter 
enables to recover, at the TD-MCSCF level, excitation energies which are rather close to 
the FCI ones as shown in Fig. [2] (a). Near the equilibrium distance the excited 
states (n = 2, . . . , 5), as described by the TD-MCSCF model, are dominated by the singly 
excited configuration (lag) 1 (nag) 1 . This statement holds at the TD-MC-srLDA level for 
the first three lowest 1 S+ excited states, which explains why, in the light of Sec. 



IV A 



the 



corresponding excitation energies are very close to the TD-HF-srLDA ones, as shown in Fig. [2] 
(b). On the other hand, the TD-MC-srLDA 5*1]+ state is dominated by the doubly excited 
configuration (la u ) 2 already at the equilibrium distance. This is illustrated by the large 
difference between the TD-HF-srLDA and TD-MC-srLDA excitation energies at 1.4 a.u. 
Note that, at the TD-MCSCF level, the doubly excited character of the 5 X E+ state appears 
only at 2 a.u., when the slope of the potential curve suddenly changes. Upon bond stretching, 
the avoided crossing of the TD-MC-srLDA 5 X S+ and 4 1 S+ states occurs around 2 a.u. The 
latter gains a doubly excited character which explains the deviation of its excitation energy 
from the TD-HF-srLDA one. Note that regular TD-LDA as well as TD-HF-srLDA methods 
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completely miss the 5 X S+ state for stretched geometries. Returning to the 4 X E+ state, the 
TD-MC-srLDA exhibits an avoided crossing with the 3 X E+ state around 2.5 a.u. The latter 
is then dominated by the doubly excited configuration (lcr u ) 2 . As a consequence the TD- 
HF-srLDA and TD-MC-srLDA 3 X S^ excitation energies differ significantly when stretching 
the bond beyond 2.5 a.u. Another avoided crossing occurs, at 3 a.u. in the TD-MC-srLDA 
model, between the 3 X E+ and 2 1 £+ states. While regular TD-LDA as well as TD-HF-srLDA 
keep on describing the lowest excited state as singly excited for large bond distances, 
the TD-MC-srLDA method is able to catch the double excitation (lcr 9 ) 2 — » (1<j u ) 2 , which 
is included in the active space. For analysis purposes, TD-MC-srLDA excitation energies 
have been computed when /i = (see Fig. [3]). As shown in the Appendix, even though 
the ground-state MC-srLDA wave function reduces to the regular KS-LDA determinant, 
the TD-MC-srLDA (/i = 0) method recovers not only the TD-LDA spectrum but also the 
double excitation (ler 9 ) 2 — > (la u ) 2 whose excitation energy is found to be twice the KS 
orbital energy difference 2(e au — e a ). This explains the crossings in Fig. |3| As pointed 



out in Sec. |IID[ the double excitation should in principle be disregarded since it does not 
correspond to a pole of the density. Still, it is interesting to note that the TD-MC-srLDA 
(ft = 0.4) avoided crossings shown in Fig. [2] (b) originate from this unphysical solution. Let 
us stress that, for fi = 0.4, the double excitation described within TD-MC-srLDA is, on the 
other hand, physical. This is ensured by non-zero orbital rotation coefficients (not shown) 
in the corresponding linear response vector as well as the multi-configurational character of 



the ground-state MC-srLDA wave function in stretched geometries. Finally, as in Sec. |IV A 



we again conclude that the inaccuracy of the TD-MC-srLDA (/i=0.4) method is essentially 
due to the approximate short-range potential and kernel used. 



C. Singlet excited states of the beryllium atom 

Singlet excitation energies of the beryllium atom have been computed at the TD-MC- 
srDFT level, using the srPBEgws functional and the 2s2p active space, when varying the 
/i parameter. Results are shown in Figs. [4] (a) and [4] (b). Comparison is made with the 
linear response TD-CCSD results which are used as reference. Let us first mention that, 
in the fi = limit of the TD-MC-srPBEgws model, standard TD-PBE excitation energies 
are recovered. Interestingly the l l D state, which corresponds to the double excitation 
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(2s) 2 — > (2p) 2 and which is absent in the standard TD-PBE spectrum, is also described 
at the TD-MC-srPBEgws level even when fi = 0. In this particular case, the excitation 
energy actually equals twice the KS-PBE orbital energy difference 2(e 2p — ^2s)- As shown 
in the Appendix, this is due to the fact that the double excitation is included in the active 
space and can, therefore, be treated explicitly. In this case, the excitation energy turns 
out to be equal the TD-CCSD one. Nevertheless, when /i is strictly equal to zero, this 



double excitation is unphysical as mentioned in Sec. |II D| For non-zero fi values, even small 
ones, this is not the case anymore since the ground-state MC-srPBEgws wave function is 
multi-configurational and orbital rotations appear in the linear response vector. Let us 
mention that, at the TD-MC-srLDA level, we obtained a l^D excitation energy of 7.17 
eV when fi = 0.1 which corresponds to an absolute error of 0.02 eV. In the TD-DMFT- 
srLDA (fi = 0.1) method of Pernal [23], the absolute error is much larger (2 eV). (We note 
that Pernal used a different srLDA functional than we used, namely the one by Paziani 
et al. [BD], but we do not expect this to be important for the absolute error.) As pointed 
out by the author this error might be due to the approximate density-matrix functional 
used for describing the long-range interaction while, in the TD-MC-srDFT model, the latter 
is explicitly treated at the MCSCF level. When increasing the fi parameter, the TD-MC- 
srPBEgws \ X D excitation energy deviates slightly from the TD-CCSD one, by less than 0.15 
eV. Note that the TD-MC-srLDA \ X D excitation energy was found to be equal to 6.96 eV 
for fi = 0.4 a.u., which corresponds to an absolute deviation of 0.23 eV. The absolute error 
obtained at the TD-DMFT-srLDA (/i = 0.4) level by Pernal [23] is much larger (3.38 eV), 
which we again attribute to the approximate density-matrix functional used for modeling the 
long-range interaction. The other excitation energies are all underestimated at the regular 
TD-PBE level. Those can be improved significantly when fi becomes larger than 0.2 a.u. 
The difference between the TD-HF-srPBEgws and TD-MC-srPBEgws excitation energies 
increases then with fi since electron correlation is transferred from the short-range DFT 
part to the long-range MCSCF part of the linear response equations. We note that, for 
most of the states, TD-MC-srPBEgws can perform better than TD-HF-srPBEgws and the 
regular TD-MCSCF model, which is recovered in the fi — > +oo limit, but different fi values 
should then be used. As extensively discussed in Ref. [21] it is more appealing, in order to 
have a general method, to set the fi parameter to a fixed value, using prescriptions for an 
optimal treatment of correlation for example, and improve the accuracy of the short-range 
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functionals. Work is in progress in this direction. 



D. ferrocene 

A major objective stimulating the development of the TD-MC-srDFT model is to pave 
the way for an efficient computational tool of predictive power that allows us to study 
excitation spectra of large, complex molecular systems where for example local and charge- 
transfer (CT) excitations often co-occur. As an illustrative example we computed the low- 
lying singlet excitation spectrum of ferrocene for which experimental solution data [6TJ [62] 
as well as ample theoretical ab initio gas-phase reference values [121 EH I63ti65] are available. 

A ground-state HF calculation for ferrocene with the Cp-rings in D 5 ^ symmetric con- 
figuration around the d 6 Fe(II) center yields the following energetic order of iron-centered 
high-lying occupied (non-bonding) and anti-bonding LUMO <i-orbitals: 

o^((i z 2, d a ) ~ €2^dx 2 —y 2 i d X yi d$) <C C]_ (d xz , dy Z , d n ). 

In addition, we find as highest occupied MOs (doubly degenerate HOMO orbitals) predom- 
inantly Cp-centered, bonding 7r-orbitals. In summary, the valence configuration reads as: 
(«'i) 2 (4) 4 (el) 4 (el?. 

Based on this ground state occupation we can expect a manifold of three doubly- 
degenerate (E';, E£, E'O singlet metal-centered valence transitions which will be examined in 



Section IV D 1 as well as ligand-to-metal CT excitations in the ultra-violet visible spectrum 



to be discussed in Section IV D 2 In particular in the latter case it will be interesting to 
examine whether our new TD-MC-srDFT scheme can improve TD-DFT/CAM-B3LYP. In 
TD-DFT/CAM-B3LYP both long- and short-range correlations are treated in DFT whereas 
in TD-MC-srDFT the long-range correlation is handled with MCSCF. 

All experimental data available to us is recorded in solution, and a fair comparison of the 
computed data would thus require to account for solvent effects which was beyond the scope 
of the present investigation. Work in that direction is currently in progress |66j. Hence, we 
assume that the metal-centered transitions are rather insensitive to solvent effects due to 
their local character whereas the CT transitions are more likely to be sensitive to solvation 
effects because of potentially strong ligand-solvent interactions. 

Before embarking on details of the excited state spectrum compiled in Table |T] we com- 
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mence with a brief examination of the geometry dependence of our data. For this purpose, 
we performed TD-DFT/CAM-B3LYP calculations with ferrocene (i) in a staggered con- 
formation with D 5 d symmetry (for geometrical data see Table S2 in the Supplementary 
Material [15]) and (ii) in a second eclipsed conformation with geometrical parameters taken 
from experiment [671 - 169] which was used by Ishimura et al. [19] in their symmetry adapted 
cluster-configuration interaction (SACCI) excited state calculations (see Table S3 in the Sup- 
plementary Material [15])- We find that the lowest-lying singlet excited states are hardly 
affected by adapting to a staggered conformation. The d-d excitation energies are lowered to 
a minor degree to 2.29 eV, 2.62 eV and 3.43 eV with a maximum deviation of A max = -0.04 
eV from the reference TD-DFT/CAM-B3LYP values (Ds/j-symmetric configuration) listed 
in Table [T] Larger deviations with A max = -0.10 eV are observed for the second eclipsed 
configuration of ferrocene where the d-d excitation energies read as 2.23 eV, 2.56 eV and 
3.37 eV. The latter should be taken into account when comparing our final results with the 
SACCI data given in Table [I] 

1. Metal- centered d-d transitions 

Table |T] comprises in the upper part the spin-allowed singlet d-d transition energies. For 
further analysis we summarize in Table [TT] the relative energy gaps Ai and A 2 , respectively 
which are determined from the excitation energy differences between E 2 ' and the lower 
state and the higher-lying singlet d-d transition E'/ and Eg. 

Turning first to the transition energies computed at the TD-HF level of theory reveals 
strikingly the importance of accounting for electron correlation. We find deviations from 
experiment by more than 1.8 eV for the lower E x ' and E 2 ' states. The comparatively close 
match of the second E^' transition of 3.39 eV with the experimental reference of 3.82 eV 
is fortuitous. The TD-MCSCF approach in contrast allows to recover a significant part of 
the dynamic electron correlation in the excited states since all <i-orbitals involved in the 
low-lying d-d transitions are included in the active CAS(10,10) space. This indeed leads to 
overall, greatly improved absolute excitation energies. The third excited state E'/, however, 
is located at only 3.15 eV giving rise to a too close gap A 2 = 0.19 eV and in addition, Ai 
hardly changes compared to TD-HF. 

The excited state data derived from regular TD-DFT calculations listed in Tables [T] and 
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\U\ follow general tendencies that can be summarized as follows. Introducing a fraction of 
full-range exact HF exchange like in B3LYP (20 %) and CAM-B3LYP (19 %) compared to 
the LDA and PBE functionals correlates with decreasing excitation energies in contrast to 
increasing energetic gaps Aj^. In general, this results in a closer agreement with experiment 
using hybrid functionals. Hence, we may order the accuracy of the functionals studied here as 
CAM-B3LYP ps B3LYP > PBE > LDA. As one could expect for the rather spatially localized 
d-d transitions the hybrid B3LYP functional performs similarly well compared to the more 
sophisticated CAM-B3LYP functional which includes an additional 46 % of long-range HF 
exchange. Let us now consider the excitation energies and their energetic separations A x 
and A2 computed by means of TD-HF-srDFT. Using 100% of long-range HF exchange, 
irrespective of the choice of srLDA or srPBE functionals, reduces the transition energies 
and enlarges their relative gaps Ai and A2 in comparison to the regular TD-DFT/LDA 
and TD-DFT/PBE data, respectively. This is not only in line with the findings for (CAM- 
)B3LYP relative to LDA and PBE but results also in a good agreement of TD-HF-srDFT 
with the hybrid functionals. Clearly, \i = 0.4 a.u. is large enough to introduce a non- 
negligible part of HF exchange, since long-range effects do not appear to be significant in 
the d-d excitations. The fact that TD-MC-srPBEgws and TD-HF-srPBEgws results are 
not equal elucidates that a non-negligible part of the two-electron interaction is indeed 
assigned to the long-range interaction and thus treated by the MCSCF approach within the 
CAS(10,10) active space. In accordance with Fromager et al. [21], total energy differences 
on the order of 10~ 2 a.u. between the ground-state HF-srPBEgws and MC-srPBEgws wave 
functions, respectively, corroborate our conclusion. In passing, we note that our TD-MC- 
srPBEgws data are in very good agreement with experiment with respect to both excitation 
energies and the relative energetic spacings Ai^ of the three lowest electronic singlet d-d 
transitions. The largest deviation to experiment is for the third excited state (located at 3.82 
eV) where all approaches studied here yield too low transition energies except for SACCI 
which, not taking into account further geometry effects {vide supra), slightly overshoots by 
+0.21 eV. In summary, based on the experimental references, our new TD-MC-srPBEgws 
method shows significantly enhanced performance compared to the multi-reference wave 
function SACCI approach which yields only a moderate agreement with experiment both 
with respect to transition energies and their energetic separation A 12 . Moreover, TD-MC- 
srPBEgws provides a valuable alternative to the otherwise well-performing hybrid (CAM- 
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)B3LYP functionals as well as the model potential LB94 which was explored in the recent 
work by Scuppa and co-workers [50] . 

2. Ligand —> Metal charge-transfer excitations 

We report in the lower part of Table [I] CT excitation energies computed at the TD-HF, 
TD-DFT level employing LDA, GGA, B3LYP and CAM-B3LYP functionals along with our 
TD-HF/MC-srDFT results. The experimental data refers to the maximum peak positions 
in the measured absorption spectrum of ferrocene. It is readily seen that the TD-HF CT 
excitation energies are unacceptably high whereas the TD-DFT/LDA and TD-DFT/PBE 
transitions are too close to the lower end of the visible spectrum. With the hybrid B3LYP 
functional, on the other hand, deviations from experiment are reduced to -0.41 and -0.48 eV 
for the first and second CT transitions, respectively, yielding excitation energies of 5.41 eV 
and 5.72 eV. This clearly indicates the importance of using an exact exchange energy, or at 
least a fraction of it, with regard to a more general description of CT states where the extent 
of exact exchange determines the partial correction of the "long-range" self-interaction error 
intrinsically present for pure LDA and GGA-type functionals. The appropriate combination 
with a suitable correlation functional then leads to satisfying results as seen here for B3LYP. 
Bearing this in mind, the superior performance of CAM-B3LYP compared to B3LYP with 
respect to the ligand-to-metal CT excitations might not be unexpected as CAM-B3LYP 
comprises in addition to the 19% of full-range HF exchange energy (B3LYP: 20%), 46% of 
the long-range exact exchange energy. Although only a fraction of long-range HF exchange is 
used, fortuitous error cancellations in combination with an approximate correlation density- 
functional may lead to excellent results, where CAM-B3LYP excitation energies at the triple- 
£ basis set level nearly coincide with the experimental solution data. Turning to the oscillator 
strengths, they are too weak in all our calculations, which is probably a combination of basis 
set effects and the solvent effects on the experiments. 

A closer inspection of the CT transitions in Table [T] further reveals that both TD-DFT/LDA 
and TD-HF predict the CT transitions in reversed order. This feature remains also in the 
picture predicted by our TD-MCSCF response calculations where the description of long- 
range and short-range dynamic correlation within the CAS(10,10) space vastly improves the 
excitation energies, lowering them to 5.50 eV for the Ag and 6.24 eV for the E' x state while, 
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at large, their relative gap (in reversed order) enlarges to 0.74 eV compared to 0.38 eV in 
experiment. The extended active orbital space treated within the SACCI approach, on the 
other hand, seems to capture most of the important contributions to the CT transitions, thus 
illustrating the importance of short-range dynamic correlation. Peaks are predicted closely 
spaced at 6.34 eV and 6.43 eV, respectively, though. Similar conclusions can be drawn for the 
TD-HF-srDFT model with either the srLDA or srPBEgws functional. The straightforward 
inclusion of short-range dynamic correlation by means of the srDFT functional, using \x = 0.4 
a.u., suffices in this context to recover the correct CT state order from the /i — > +oo pure 



wave-function-based TD-HF limit. In contrast to the local d-d transitions (see Sec. IV D 1) 
we observe for the first CT state a large difference between the TD-HF-srDFT and TD- 
DFT/CAM-B3LYP approaches which we ascribe mainly to different fractions of long-range 
HF exchange, which amounts to only 65% in CAM-B3LYP, and more importantly, to missing 
long-range correlation in TD-HF-srDFT. The impact of describing the latter correlation type 
explicitly, within the MCSCF model, is highlighted by our TD-MC-srPBEgws results. It 
lowers the first E' x transition by 0.35 eV compared to TD-HF-srPBEgws while the second 
main CT peak is pushed upwards to 7.04 eV, thus overshooting the experimental reference 
by « 0.8 eV. A detailed analysis shows that the second CT transition is, within both the 
TD-MCSCF and TD-MC-srPBEgws approaches, fully described by orbital rotations whereas 
the optically allowed E x transition is largely obtained by configuration rotations. However, 
one cannot make any final conclusions about the quality of the TD-MC-srPBEgws results 
compared to experiments before the solvent effects have been included in the calculations. 
Work on this is currently in progress. Computing high-level gas phase CC excitation energies 
would also be interesting for benchmarking. 

To sum up, our new proposed TD-MC-srDFT model has promising potential to accurately 
describe not only local excitations as shown for example for the d-d transitions in Sec. |IVD 1 



but also to simultaneously model charge transfer excitations and their corresponding relative 
oscillator strengths. It exhibits in particular improved performance compared to both TD- 
DFT based on pure density-functionals and TD-MCSCF. In contrast to the latter approach, 
the srPBEgws functional combined with a long-range TD-MCSCF treatment is able to 
describe short-range dynamic correlations and enables to get the two CT states in predicted 
order. The relative energy difference between the two CT states seems too large, though, 
when compared to experimental data. This might be due to self-interactions errors induced 
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by the employed short-range functionals where the formers could be reduced choosing a 
different separation of short-range exchange and correlation energies. Work is currently 
in progress in this direction. Moreover, the inclusion of solvation effects by means of a 
polarizable embedded framework into the TD-MC-srDFT and TD-DFT response approaches 
[66] of the DALT0N2011 program framework will facilitate future comparison to experimental 
works recorded in solution. Finally, we note that the TD-MC-srDFT scheme seems to 
be a viable alternative to TD-DFT/CAM-B3LYP where, in the former case, long-range 
correlation effects are described at the MCSCF level. Further calibration studies focusing 
particularly on such a direct comparison are being prepared. 



V. CONCLUSIONS 



The extension of multi-determinant range-separated DFT to the time-dependent regime 
has been investigated. Following Vignale [251 126] . an exact variational formulation was 
obtained. Various approximate schemes can then be formulated, depending on the choice of 
the post-HF method that is used for describing the long-range interaction, and the choice 
of the exchange-correlation functional used for treating the short-range interaction. In this 
work, the combination of long-range MCSCF with short-range adiabatic LDA and GGA 
was considered. The corresponding linear response scheme, referred to as TD-MC-srDFT, 
was then derived within Floquet theory. Numerical results obtained for the singlet excited 
states of the stretched H2 molecule and Be show that TD-MC-srDFT can, in contrast to 
regular TD-DFT, describe double excitations even though the adiabatic approximation is 
used for the short-range exchange-correlation density-functional contributions. This is made 
possible by the MCSCF-type parametrization of the wave function where, in addition to 
orbital rotations, configuration rotations inside the active space can be treated explicitly. 
For the commonly used \x = 0.4 range-separation parameter value, excitation energies are, 
in most cases, much better described with TD-MC-srDFT than regular TD-DFT. Huge 
improvements are also observed, when comparing with the TD-DMFT-srLDA results of 
Pernal [23], in the description of the \ X D doubly-excited state of Be and the state 
of the stretched H 2 molecule. Still, in the latter case, excitation energies are significantly 
underestimated. The error is due to the approximate short-range LDA and GGA potential 
and kernel that are used. Increasing \x improves on the accuracy but then some part of the 
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short-range electron correlation is transferred from DFT to the MCSCF. This is of course not 
appealing, in terms of computational cost, for larger scale calculations. As an alternative, a 
different decomposition of the short-range exchange and correlation energies [SHHSO] could 
be used while keeping the fi parameter set to 0.4 a.u. Work is currently in progress in 
this direction. The accuracy of the TD-MC-srDFT approach using a short-range GGA 
functional was further examined for the low-lying singlet excited states of the d G metallocene 
ferrocene. Excitation energies as well as energetic state separations are in overall good 
agreement with experiment for both valence and charge-transfer excitations apart from a 
single outlier. TD-MC-srDFT generally outperforms traditional wave-function and TD-DFT 
approaches with exception of TD-DFT/CAM-B3LYP that yields comparable results with 
a better performance for the charge-transfer states. The latter functional was specifically 
designed to remedy known shortcomings of TD-DFT [18] for this excitation class, though. 
Extensive works with the overall aim to shed further light on weaknesses and strengths of 
our proposed TD-MC-srDFT approach with respect to TD-DFT/CAM-B3LYP are therefore 
in progress. Particular focus is in this context laid on model excited-state charge-transfer 
compounds such as peptides, carotenoides [70] and transition metal complexes [71] . 
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APPENDIX: TWO-STATE TD-MC-srDFT MODEL IN THE /x = LIMIT 

Let us consider a minimal active space consisting of two Slater determinants \a 2 ) and 
\b 2 ) representing doubly-occupied <fr a and (pt orbitals, respectively. The former determinant 
corresponds to the ground-state KS determinant which is recovered at the MC-srDFT level 



in the /i = limit (see Sec. II A). For simplicity, orbital rotations are not considered in the 
following. According to Eqs. (7) and (8) in Ref. |41| . the long-range Hessian and the 
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metric S^ 2 ^ introduced in Eqs. (38) and (39) become, when \x = 



where 



A 



B c 



A 



E 



[2]0 













A 








A 




B°* 


A * 




-A * 





(a 2 \[R a ,, [H KS , Rl]]\a 2 ) (a 2 \[R a2 , [tf KS , £[ 2 ]]|a 2 > 

(a 2 |[i^if K ^ 2 ]]|a 2 > (a 2 |[^,[^ KS ,^ 2 ]]|a 2 )_ 

(a 2 !^, [tf KS , R a ,]]\a 2 ) (a 2 \[R a2 , [H KS , R b ,]]\a 2 ) 

(a 2 \[R b *,[H KS ,R a2 ))\a 2 ) (a 2 | [jfc, [H KS , R b *)) |a 2 ) 
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(a 2 \[R b 2 } R a 2]\a 2 ) (a 2 \ R b *\ \a 2 ) 



(Al) 



(A2) 



Rl = \a 2 )(a 2 \, R\ 2 = \b 2 )(a 2 \, 
and H KS = T+1^ e +VH XC [n a 2] denotes the non-interacting KS Hamiltonian. Since H KS \a 2 ) 



2e a \a 2 ) and H KS \b 2 ) = 2e b \b 2 ), the matrices in Eq. (A2j) can be simplified as follows 



2(e b -e a ) 



A 



B° = A = 0, 



(A3) 




1 



In addition, the gradient density vector defined in Eq. ( |37[ ) becomes in the /i = limit: 



(b 2 \h(r)\a 2 ) 


-(a 2 |n(r)|6 2 ) 



n 



W°(r) 



0. 



(A4) 
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since \a 2 ) and \b 2 ) differ by a double excitation. As a result the kernel contribution to the 



Hessian defined in Eq. (40) equals zero when /i = so that the total Hessian is simply Eq 



[2]0 



We thus conclude from Eqs. (Al) and (A3) that the double excitation a 2 —> b 2 is described 
within the TD-MC-srDFT model in the /i = limit and the corresponding excitation energy 
is twice the KS orbital energy difference 2{e\ ) — e a ). In this respect, the TD-MC-srDFT model 
does not reduce to standard TD-DFT, as the TD-HF-srDFT model does, in the /i = limit. 
Still, in this example, the double excitation is a pole of the TD-MC-srDFT wave function 
but not of the density since the linear response contribution to the density n^°^(r)X(uj) 



always equals zero, according to Eq. (A4). From this point of view, the double excitation 
obtained at the TD-MC-srDFT (/i = 0) level is unphysical. 
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FIGURE CAPTIONS 



Figure [TJ First (black), second (red), third (green) and fourth (blue) excitation en- 
ergies in H 2 along the bond breaking coordinate calculated with (a) standard TD- 
MCSCF and TD-LDA methods (b) TD-HF-srLDA and TD-MC-srLDA schemes. Com- 
parison is made with FCI results. Each type of line corresponds to a given method. 
The \i parameter was set to 0.4 a.u. The minimal active space \a g \o u was used. Basis 
set is aug-cc-pVQZ. 

Figure [2} First (black), second (red), third (green) and fourth (blue) X E+ excitation en- 
ergies in H2 along the bond breaking coordinate calculated with (a) standard TD- 
MCSCF and TD-LDA methods (b) TD-HF-srLDA and TD-MC-srLDA schemes. Com- 
parison is made with FCI results. Each type of line corresponds to a given method. 
The \i parameter was set to 0.4 a.u. The minimal active space la g la u was used. Basis 
set is aug-cc-pVQZ. 

Figure [3f First (black), second (red), third (green) and fourth (blue) 1 S+ excitation en- 
ergies in H 2 along the bond breaking coordinate obtained with the TD-MC-srLDA 
method, setting fi = 0, and compared to standard TD-LDA results. For analysis pur- 
poses, twice the la u and la g KS-LDA orbital energy difference is also plotted. Each 
type of line or point corresponds to a given method. The minimal active space \a g \a u 
was used. Basis set is aug-cc-pVQZ. 

Figure [4f Excitation energies obtained in Be when varying fi at the TD-HF-srPBEgws 
and TD-MC-srPBEgws levels for (a) the 1 X P (black), 2 1 S (red), l l D (blue) and 2 1 D 
(purple) singlet states (b) the 2 1 P (grey), 3 1 P (cyan) and (green) singlet states. 
Comparison is made with TD-CCSD results. Each type of line corresponds to a given 
method. Twice the KS-PBE orbital energy difference 2(e 2p — £2s) is also plotted for 
analysis purposes. The minimal active space 2s2p was used. Basis set is aug-cc-pVQZ. 
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FIG. 1: Fromager et al, Journal of Chemical Physics 
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FIG. 2: Fromager et al, Journal of Chemical Physics 
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FIG. 3: Fromager et al, Journal of Chemical Physics 
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FIG. 4: Fromager et al, Journal of Chemical Physics 
(a) Be [1 1 P, 2 1 S, 1 1 D, 2 1 D] 
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TABLE CAPTIONS 



Table [H Key spin-allowed singlet Fe d-d and lowest optically active singlet transition en- 
ergies (eV) in ferrocene. Oscillator strength / are reported in the row below the 
charge-transfer state excitation energies for selected, exemplary cases. The abbrevia- 
tion "TD-" has been omitted for compactness in the table header. 

Table [TTJ Energy gaps (in eV) between the lowest three spin-allowed singlet Fe d-d tran- 
sition energies in ferrocene obtained from TD-MC-srDFT calculations, and compared 
with other ab initio and experimental data. Aj refers to the energetic difference 
between the lower E x and Ej states while A2 quantifies the difference in transition 
energies between the upper E'/ and E^ states, respectively. The abbreviation "TD-" 
has been omitted for compactness in the table header. 



39 



fa 



fa 



< 
u 



a 
x 
fa 



fa 

m 
fa 

!h 

;/3 

O 
fa 

fa 

<J 

Q 
fa 



fa 

$ 
co 
fa 
fa 

% 

CO 

fa 

S 1 
OS 

fa 
fa 

fa 
fa 
fa 

< 

Q 
fa 



O 
O 
<J 
en 

fa 

o 

m 
O 



fa 



CJ 

to 
T3 



o 


00 


CM 






•y_ 


CM 




CO 


02 


'00 


O 




'00 




cn 


cm 


CO 


CM 




! — 1 




t> 




cm 


cm 


CO 










l>- 


CM 


cm 


cm 


CO 


CO 


10 


10 


CO 


CD 




cm 


cm 


CO 


CO 





00 










cm 


CO 




, — 1 




00 








cm 


CO 





CO 





05 







CM 


CO 


CO 


-t 




CM 








o-i 


cm 


CO 




t- 


CO 


1— 1 


cn 





ci 


cm 


-p 


CM 





LO 


00 


05 


I — I 


cm 


cm 


CO 


CM 











CO 





1 — 1 


CO 











CM 




O 




OO 




CM 




LO 




CO 





© 



CJ 

to 
s- 

a? 

s 
£ 

u 



05 

OO 



CM 



■OO 
CN 



LO 

'00 



LO 



o 



CO 



CM 



G5 



»3 



2- 



o 



CO 



CO 
CO 



LO 

CM 
CO 



CM 



CO 



00 



CO 



o 

LO 



O? 



o 

CO 



fa 



fa fa 



fa 



40 



fa 



fa 



< 
u 



Pi 

cu 



fa 



fa 
pq 
fa 



fa 
fa 
fa 

< 
Q 
fa 



fa 

CO 

fa 
fa 

CO 

fa 

-p 

fa 
fa 
fa 
fa 
fa 

< 
fa 
fa 



o 

O 
< 

CO 

fa 

O 

CO 

O 



fa' 



>> 
aO 
t-i 
o 

fa 



QC 






CM 


00 




o 


o 






CM 




CM 


CO 




o 


o 






CM 




CO 


CD 




o 


o 






LO 




CM 


LO 










CM 


o 




CO 


00 




o 


o 




£^ 


00 




CM 














CO 




, — 1 


LO 




o 


o 




CO 






1 — 1 


LO 




o 


o 


1— 1 






CM 


LO 


CO 


|co| 


o 


LO 




o 


o 


* l ' 

cd 














cd 


CD 


of 


1 — 1 


l> 


fe 


o 


1 — 1 


bD 






fa 






fa 






fa 




CI 




1 — 1 


1 — 1 


CO 


o 


o 













o 


|lo| 


I — 1 


CO 




o 


CM 


CD 






« 












jmj 








< 


< 


CD 












a 



41 



